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ABSTRACT 


The  topic  of  this  thesis  is  the  flow  of  pseudo-plastic 
materials.  A  discussion  of  the  empirical  model  used  to 
describe  the  behavior  of  this  type  of  material  is  presented, 
and  objections  to  its  use  are  discussed.  It  is  shown  that 
the  incompressible  Newtonian  fluid  and  the  rigid-perf ectly 
plastic  von  Mises  solid  are  special  cases  of  this  model. 

Two  problems,  plane  flow  in  a  converging  channel  and  axi- 
symmetric  flow  in  a  cone,  are  solved  for  various  "degrees" 
of  non-Newtonian  behavior.  Body  forces  and  thermal  effects 
are  neglected,  and  the  flow  is  assumed  to  be  fully  developed. 
A  numerical  technique  developed  by  Clutter  and  Smith  to  solve 
boundary-layer  problems  is  adapted  to  solve  the  non-linear 
ordinary  differential  equations  that  govern  the  flow  problems 

considered  here. 

To  obtain  the  solution  for  fully  developed  flow,  the 
streamlines  are  assumed  to  be  radial,  meeting  at  the  virtual 
apex  of  the  converging  channel  or  cone.  An  attempt  was  made, 
with  no  success,  to  consider  entrance  effects  for  the  flow 
of  a  Newtonian  fluid  in  a  cone;  this,  and  the  problem  of 
entrance  effects  for  the  pseudo-plastic  material,  is  dis¬ 


cussed  . 


. 
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NOTATION 


Where  cartesian  tensor  notation  is  used  the  usual 
summation  convention  is  employed 


o  .  . 

1J 

cartesian  components  of  the  stress  tensor 

referred  to  axes  OX. 

1 

o  .  .  , 
ij  ,k 

comma  denotes  partial  differentiation  with 
respect  to  x^ 

v  . 

1 

components  of  the  velocity  vector 

d  .  . 

components  of  the  rate  of  strain  tensor 

s  .  . 

components  of  the  deviatoric  stress  tensor 

V 

a  rheological  parameter  in  the  power  law  model 
and  the  viscosity  of  a  Newtonian  fluid 

h 

d..,  a  strain  rate  invariant 
i  l 

4  d  d..,  a  strain  rate  invariant 

2  ij  13 

h 

4  d  .  d  d  a  strain  rate  invariant 

3  im  mn  ni 

k 

shear  yield  stress  for  von  Mises  perfectly 
plastic  solid 

a 

semi-angle  of  the  channel  or  cone 

-P 

the  hydrostatic  part  of  the  stress  tensor 

Q 

volume  flow  rate 

P 

the  hydrostatic  stress 

a!’a2 

rheological  parameters  of  the  Stokesian  fluid, 
called  the  apparent  viscosity  and  cross  viscosity, 
or  viscous  normal  stress  coefficient,  respectively 

' 


work  function 


Stokes ' 


stream  function 


CHAPTER  I 


INTRODUCTION 

Rheology  is  a  branch  of  continuum  mechanics  concerned 
with  the  study  of  anything  that  flows.  "A  continuum  is  a 
(hypothetical)  structureless  substance  to  which,  at  any 
point,  we  can  assign  kinematical  or  dynamical  variables 
which  are  continuous  functions  of  the  spatial  coordinates 
of  that  point"  (1) .  The  solution  of  a  flow  problem  must, 
therefore,  begin  with  a  consideration  of  the  basic  equations 
of  continuum  mechanics. 

Continuum  mechanics  is  based  on  the  momentum  equations, 
the  equations  of  continuity  and  energy,  and  constitutive 
equations.  Assuming  there  are  no  body  moments,  the  stress 
tensor  is  symmetric  with  six  independent  components.  The 
energy  equation  is  not  considered  since  thermal  effects  are 
neglected  in  this  thesis.  Neglecting  body  forces  and  con¬ 
sidering  incompressible  quasi-static  flow,  the  equations  of 
motion  and  continuity  are,  respectively, 


>  j 


0  , 


(1.1) 


and 


v 


i 


(1.2) 


Equations  (1.1)  and  (1.2)  provide  four  equations  for  the 
nine  unknowns  and  v ^ . 


This  is  all  that  can  be  determined 


2 


from  kinematics  and  dynamics.  Constitutive  equations  des¬ 
cribe  the  intrinsic  behavior  of  a  material  in  terms  of 
stress,  stress  rates,  deformation,  and  deformation  rates, 
and  provide  the  additional  equations  needed  to  solve  problems. 
The  requirements  for  constitutive  equations  are  discussed 
by  Aris  (2).  The  constitutive  equations  used  in  this  thesis 
are  discussed  in  the  next  chapter.  Equations  (1.1),  (1.2), 

and  the  constitutive  equations,  along  with  the  appropriate 
boundary  conditions,  are  sufficient  to  solve  the  problems 
considered  in  this  thesis. 


FIGURE  1.1 

STRESS-STRAIN  RELATIONSHIPS 
IN  SIMPLE  SHEAR 


3 


1  .  1  Pseudo-Plastic,  ity 

Fluids  which  have  an  apparent  viscosity  which  increases 
with  increasing  strain  rate,  as  measured  in  simple  shear 
flow,  are  called  dilatant,  and  those  with  an  apparent  vis¬ 
cosity  which  decreases  with  increasing  strain  rate  are 
called  pseudo-plastic .  The  name  dilatant  is  actually  a 
misnomer,  since  it  implies  an  increase  in  volume.  Brodkey 
(3)  reserves  this  terminology  for  materials  which  actually 
dilate,  and  calls  the  type  of  behavior  considered  here 
shear-thickening.  A  qualitative  picture  of  the  behavior 
of  these  two  types  of  fluids  is  given  in  Figure  1.1. 

The  model  used  in  this  thesis  to  describe  pseudo¬ 
plastic  flow  is  the  power  law  model, 


n  - 1 

S±j  =  2y  ( 2 1 2 )  2  d±j  , 


(1.3) 


where  the  d-^j  are  the  components  of  the  strain  rate  tensor, 


di j  2  i >  j  +  V J  * i ^  ’ 


and  I2  is  a  strain  rate  invariant, 


This  model  is  a  generalization  of  the  model  suggested  by 
0  s  t  w  a  1 d  for  simple  shear  flow,  and  represents  dilatant  t  L  u i d  s 
for  n  >  1,  and  pseudo -plastic  fluids  for  n  <  1.  For  n  =  1 


-fc 

' 
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it  represents  the  incompressible  Newtonian  fluid,  and  for 
n  =  0  the  r igid-per f ect ly  plastic  von  Mises  solid,  with 
constitutive  equation 


(1.4) 


where 


X  = 


and 

k  =  /2  y 

Only  pseudo-plastic  materials,  and  the  special  cases  of  the 
Newtonian  fluid  and  the  von  Mises  solid,  are  considered  in 

this  thesis . 

1 . 2  Known  Solutions 

The  problems  considered  are  plane  flow  in  a  converging 
channel  and  axi-symmetric  flow  in  a  cone.  Solutions  to 
these  two  problems  have  been  obtained  for  the  special  cases 
of  the  Newtonian  fluid  and  the  r igid-per fectly  plastic  von 
Mises  so  lid . 

For  the  Newtonian  fluid,  the  effect  of  the  inclusion 
of  inertia  terms  has  been  considered  for  both  problems. 
Rosenhead  (4)  obtained  solutions  in  terms  of  elliptic  func¬ 
tions  for  the  plane  flow  problem  with  inertia  effects  included 
These  solutions  indicate  that,  for  a  given  semi-angle  and 
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Reynolds  number,  an  infinite  number  of  velocity  profiles 
are  possible.  These  profiles  may  or  may  not  be  symmetric 
with  respect  to  the  center  line  of  the  channel,  and  inflow 
and  outflow  can  occur  at  the  same  time.  Rosenhead  suggested 
that  stability  considerations  would  probably  exclude  many 
of  these  flows,  and  that  the  flow  pattern  obtained  in  an 
experiment  would  be  determined  by  the  pressure  conditions 
at  the  inlet  and  outlet.  Hamel  (5)  showed  that  no  radial 
flow  solution  exists  for  the  axi-symmetr ic  problem  if  inertia 
terms  are  included.  Simple  radial  flow  solutions  are  obtained 
for  both  problems  when  inertia  effects  are  neglected.  These 
solutions  are  outlined  below. 

The  solution  for  the  plane,  quasi-static  flow  of  a 
Newtonian  fluid  in  a  converging  channel  with  semi -angle  a 
is  easily  found  to  be 

v  =  2Q-  sin2_a _  [i  _  (silL6-) 2  ]  I  ,  (1.5) 

r  sin2a  -  2acos2a  sina  r 

and 

=  4yQ  ,(sin2e-3)  _  (sln290-3) 

P  sin2a  -  2acos2a  r2  r^2 

where  pq  is  the  hydrostatic  pressure  at  0q >  ,  Q  is  the 

a 

Q  =  2r  /  vr  d0  , 

0 


volume  flow  rate, 
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and  y  is  the  coefficient  of  viscosity. 

The  corresponding  solution  for  the  a x i - symme t r i c  pro¬ 
blem  was  given  by  Bond  (6).  This  solution  can  be  written 
in  the  form 


- - -  [1  -  (s±n0_j2  j  1  >  (1.6) 

2tt  ( 1-cosa)  2  (l+2cosa)  sina  r2 


and 


P  -  P 


0 


(2-3sin26 )  (2-3sin20o) 


ir(l-cosa)2  (l+2cosa) 


0 


where 

a 

Q  =  2  tt  r  2  /  vr  sin0  dO 
0 


Nadai  obtained  the  stress  field  for  the  flow  of  a  rigid- 
perfectly  plastic  von  Mises  solid  in  a  converging  channel, 
with  the  assumption  that  the  shear  yield  stress  k  is  gener¬ 
ated  at  the  sides.  The  associated  velocity  field  was  found 
by  Hill  (7).  A  shear  stress  boundary  condition  at  the  channel 
wall  was  used  in  solving  this  problem,  in  contrast  to  the 
no-slip  boundary  condition  used  in  obtaining  the  solutions 
for  the  Newtonian  fluid. 

Shield  (8)  considered  the  flow  of  a  r i gi d -p er f e c t ly 
plastic  solid  in  a  cone.  Shield  did  not  obtain  a  closed 
form  solution,  but  integrated  the  equations  numerically. 


" 
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Again,  a  shear  stress  boundary  condition  was  used. 

These  known  solutions,  which  are  special  cases  of  the 
flows  studied  in  this  thesis,  provide  a  check  on  the  results 
obtained,  and  the  solutions  for  the  Newtonian  fluid  are  used 
as  a  starting  point  in  the  numerical  procedure. 

1 . 3  Aim  of  This  Thesis 

The  aim  of  this  thesis  is  to  study  pseudo-plastic  flow 
in  a  converging  channel  and  in  a  cone,  and  to  develop  an 
efficient  numerical  technique  for  solving  the  non-linear 
two-point  boundary-value  problems  for  fully  developed  flow. 


CHAPTER  II 


THE  POWER  LAW  MODEL 

In  this  chapter,  the  power  law  is  considered  as  an 
emp ir ical  appro xima tion  for  the  constitutive  equations  of 
certain  Stokesian  fluids.  This  model,  along  with  the  equa¬ 
tions  of  motion  and  continuity,  is  used  to  describe  pseudo¬ 
plastic  fl  ow  . 

2 . 1  The  Stokesian  Fluid 

The  Stokesian  Fluid  is  an  isotropic,  homogeneous, 
inelastic  fluid  which  satisfies  the  following  relations 
(9),  (10); 


ai  j  ai  j  ^r  s  )  5 

and 

a±J  (0)  =  -  P  5±j  ,  (2.1) 

where  P  is  the  hydrostatic  stress.  The  components  of  the 
stress  tensor  can  be  written  as  the  sum  of  a  deviatoric  part 
and  a  hydrostatic  part, 


G 


ij 


(2.2) 


where  the  S,--  are  components  of  die  deviatoric  stress  tensor. 


and 
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Assuming  the  components  of  the  viscous  stress  tensor  can 
be  expressed  as  polynomials  in  the  components  of  the  rate 
of  strain  tensor,  the  constitutive  equations  for  the  in¬ 
compressible  Stokesian  fluid  may  be  written  as 


<*ij  =  "  P  6ij  +  al  dij  +  a2  dim  dmj  >  (2* 

since,  by  using  the  Cay ley-Hamilton  theorem,  all  terms  of 
higher  order  can  be  written  in  terms  of  6-^j  ,  d-jj  ,  d^m  dmj  , 
and  the  two  non-vanishing  invariants  of  the  rate  of  strain 
tensor , 


12  ~  \  dij  dij  > 

and 

^3  =  3^  dim  dmn  dni  * 

The  parameters  and  a 2  are  functions  of  l£  and  I3,  and 

must  be  such  that  the  rate  of  viscous  dissipation  is  non¬ 


ai  d^j  +  a2  d-j.m  dmj  dij 


=  2a^  I2  +  3»2  I3  d 


negative. 


■ 
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The  constitutive  equations  of  the  Stokesian  fluid  are  non 
linear,  and  its  rheological  parameters  depend  on  the  strain 
rate.  Also,  these  fluids  exhibit  normal  stress  effects  if 
a 2  f  0.  Equations  (2.3)  are  the  general  form  of  the  consti¬ 
tutive  equations  for  the  incompressible  Stokesian  fluid. 

2 . 2  The  Power  Law 

The  mathematical  difficulties  encountered  in  applying 
equation  (2.3)  are  extreme.  Also,  experimentalists  have  not, 
as  yet,  been  able  to  determine  the  dependence  of  a,]_  on  I3 , 
or  to  measure  a2  unambiguously.  Consequently  simpler  rela¬ 
tions  have  been  suggested.  Some  of  these  are  given  by 
Bird  (11) . 

For  some  simple  flow  problems  normal  stress  effects 
have  no  influence  on  the  relationship  between  shear  stress 
and  shear  rate.  If  the  relationship  between  shear  stress 
and  shear  rate  alone  is  of  interest,  it  is  reasonable  to 
assume  that  a2  is  zero.  If  the  existence  of  a  work  function 


E  =  E(dij) 


is  assumed  (12)  such  that 


3E(drs) 

sij  =  Til— 


then,  for  an  incompressible  fluid. 
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E  =  E(I2,I3)  , 


and 


3  E 

9dij 


3  E 
3 1 3 


( d ik  dk  j 


ij 


For  an  incompressible  Stokesian  fluid 


sij  =  al  dij  +  a2<dik  dkj 


6  .  ,  ) 
ij  J 


where 


3E 

3I2 


al  (  x2  » 1 3  >  » 


and 


3  E 
3 1 3 


a2  ^ 1 2  > 1 3  ^  • 


The  assumption  that  ot 2  is  zero  therefore  implies  that 


E  =  E(I2)  , 


ax  =  a1(I2)  , 


and 


Hence,  the  assumption  that  the  coefficient  a  in  L  he 


p  ower 


law , 


12 


n- 1 

ax  =  2y(2I2)  2  , 

does  not  depend  on  the  third  strain  rate  invariant  has  some 
theoretical  justification. 

Apparent  viscosity  can  be  measured  using  a  Couette  or 
capillary  viscometer.  Since  I3  is  zero  for  the  flow  in  these 
instruments,  the  dependence  of  012  on  I3  can  not  be  determined 
from  these  tests.  Tanner  (13)  proposed  the  use  of  quasi¬ 
static  flow  in  a  converging  cone,  in  which  I3  is  non-zero, 
to  check  the  validity  of  the  power  law.  If  the  work  function 
exists,  the  non-zero  value  of  I3  will  not  affect  the  flow, 
and  the  velocity  profile  and  pressure  will  agree  with  that 
predicted  by  the  power  law  if  the  constants  n  and  2 y  have 
been  chosen  to  fit  the  power  law  to  the  results  obtained 
from  viscometer  tests.  Approximate  values  for  n  and  2y  for 
several  fluids  are  given  in  (14). 

Reiner  (15)  raised  three  objections  to  the  power  law 
model.  First,  for  zero  strain  rates  and  n  less  than  one, 
the  apparent  coefficient  of  viscosity  is  infinite,  and 
second,  for  infinite  strain  rates,  it  is  zero.  Since  neither 
of  these  extremes  is  of  interest  in  practical  analysis, 
these  objections  are  not  serious.  Fredrickson  (16)  shows 
that  the  power  law  can  be  made  to  fit  experimental  data  over 
a  certain  range  of  rates  of  deformation.  Thus,  despite  the 
above  objections,  the  power  law  is  valid  if  the  rates  of 
deformation  are  suitably  restricted.  Ihe  third  objection 
is  more  serious.  The  dimensions  of  the  parameter  y  depend 


. 
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on  n,  and  consequently  the  power  law  can  not  be  a  physical 
law.  However,  if  n  is  a  constant  for  a  given  fluid,  and  y 
is  determined  from  experimental  data,  then  the  power  law 
is  a  valid  empirical  relationship. 


CHAPTER  III 


PSEUDO-PLASTIC  FLOW 

In  this  chapter  the  governing  equations  for  the  flow 
problems  under  consideration  are  derived. 

3 . 1  Plane  Flow  in  a  Converging  Channel 

The  first  problem  considered  is  plane,  quasi-static 
flow  of  an  incompressible  material  in  a  converging  channel. 
Plane  polar  coordinates,  r  and  9  (Figure  3.1),  are  used, 
and  the  flow  is  assumed  to  be  radial.  The  assumption  of 
radial  flow  is  justified  since  a  non-trivial  solution  is 
obtained . 

For  radial  flow,  the  components  of  the  rate  of  strain 
tensor  are,  in  polar  coordinates, 

3  vr 

dr  =  » 


and 


3  v 


[r9 


2r  30 


The  continuity  equation  is 


3  vr  vr 
- £  +  _£  =  0  . 


(3.1) 


3r 


r 


(3.2) 
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Equation  (3.2)  requires  that 


N 

s 


FIGURE  3.1 

FLOW  IN  A  CONVERGING  CHANNEL 


vr 


g(0) 

r 


(3.3) 


Substitution  of  equation  (3.3)  into  equations  (3.1)  gives 


dy  =  -  d, 


=  _  8(6) 


and 
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dr  e 


g*  (9) 

2r2 


(3.4) 


where  the  prime  denotes  differentiation  with  respect  to  0. 
For  radial  flow,  the  equations  of  motion  become 


r 


3are 

36 


a. 


a0 


0  , 


and 


r 


9ar0  , 

3r 


3a0 

30 


+  2  ar0 


0  . 


Substitution  of  equations  (2.2)  in  these  equations  yields 


and 


3 p  _  9sr  ,  2  .  1  9Sr0 

37  "  37“  r  r  30  * 


3_p 

30 


2Sre 


+  r 


3  sr  0 
3  r 


(3.5) 


since 


S@  -  Sr  . 

The  hydrostatic  pressure  can  be  eliminated  from  equations 
(3.5).  This  gives  the  single  partial  differential  equation 

32Sr  2  1  32sr9  _  32sr9  _  3  3  S.r  9  _  Q  _  (3.6) 

2  3r36  1  39  r  302  3r2  3r 
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Substitution  of  equations  (3.4)  into  the  power  law 

yields 


=  -  S0  =  - 


3-n 


2n 


n-1 


pgi^  g2  +  (g')2i  2 


and 


•re 


1  -n 


2n 


n-1 


yg [ 4  g2  +  (g  ’ ) 2  ]  2 


for  the  deviatoric  components  of  the  stress  tensor .  Let 


m 


<0)  -  . 


(3.7) 


where  m  is  a  non-dimensional  velocity  depending  on  0  only, 
and 

V  =  g (0) 

is  a  constant  to  be  determined  from  the  volume  flow  rate. 

With  the  substitution  (3.7)  the  stress  deviators  become 


Sr  = 


-  So  =  - 


3-n 


2n 


yVn  mG  , 


and 


’re 


1-n 
2  2 
„  2n 


pVn  m'G  , 


(3.8) 
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where 


n-1 

G  =  [4m2  +  (m' ) 2 ]  2 


Substitution  of  equations  (3.8)  in  equation  (3.6) 

yields 

(m'G)"  +  4n(l-n)  m'G  +  4(2n-l)(mG)'  =  0  .  (3.9) 

Substituting  equations  (3.8)  in  equations  (3.5)  and 
integrating  gives  the  hydrostatic  pressure 

1-n 

p  =  -2_—  pVn  [-  (1-n)  mG  -  (m'G)’]  +  A  ,  (3.10) 

v  r  2n  n  2n 

where  A  is  an  arbitrary  constant. 

For  n  =  1  equations  (3.9)  and  (3.10)  reduce  to  the 
equations  for  quasi-static  flow  of  a  Newtonian  fluid  in  a 
converging  channel.  De  Vries  (17)  showed  that  equation 
(3.9),  with  n  =  0,  reduces  to  the  equation  found  by  Nadai 
for  the  flow  of  a  r igid-per f ectly  plastic  von  Mises  solid 
in  a  converging  channel.  Equation  (3.10)  can  be  rearranged 
to  eliminate  n  from  the  denominator.  Dividing  equation  (3.9) 
by  2n ,  integrating,  and  rearranging,  yields 

0 

-  (1-n)  mG  -  —  (m'G)’  =  2mG  +  2 (1-n)  /  m'G  dB  +  D  , 

n  2n  0 


19 


where  D  is  a  constant.  The  expression  for  the  hydrostatic 
pressure  becomes,  for  n  =  0  , 


0 

p  =  /2  y[2mG  +  2  /  m'G  d 0 ]  +  E(r)  +  A  . 

0 


This  is  equivalent  to  the  expression  found  for  the  pressure 
from  Nadai's  equations,  for  with  the  substitution 

is  -  -  tan  2* 

De  Vries  showed  that  equation  (3.9)  reduces  to 

i[»'=c  sec  2\p  -  1  . 

Using  this  substitution  the  expression  for  the  pressure 
becomes 


=  /2  Vl-  cos  24,  +  2  /  cos  24>  sin  2i|.  dftj  +  E(r)  +  A> 

q  c-cos  2^ 


or,  by  performing  the  integration, 


p  =  /2  y  c  ln(c-cos  2i[»)  =  E(r)  +  A  .  (3.11) 

The  function  E(r)  can  be  evaluated  since,  from  equation 


(3.10)  , 


■ 
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(|P)  =  iLX L  [  (m  '  G  )  '  -  4mG], 

3r  n=0  r 


which,  with  the  above  substitution,  becomes 


(f£) 

3r  n=0 


=  2 


/2*y 


From  equation  (3.11), 


<u> 


n=0 


E'(r)  , 


and  therefore 


E  (r )  =  2/2  yc  In  r  +  B  . 

The  pressure  now  becomes 

p  =  /2  y[c  ln(c-cos  2\p)  +  2c  In  r]  +  F  , 

which  is  the  expression  derived  from  Nadai's  equations. 

These  results  show  that  these  two  problems  are  special  cases 
of  the  flows  considered  in  this  thesis. 

The  problem  reduces  to  the  solution  of  equation  (3.9) 
subject  to  the  boundary  conditions 


m (0)  =  1  , 


. 
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m' (0)  =  0 

from  symmetry,  and  the  no  slip  boundary  condition 

m(a)  =  0  , 

where  a  is  the  semi-angle  of  the  channel.  In  general,  this 
must  be  done  numerically. 

For  n  =  De  Vries  (17)  showed  that  equation  (3.9) 
can  be  integrated  analytically  for  m'G.  The  solution  obtained 
is 


m'G  =  c  s in0  , 

and  the  shear  stress  becomes 

1  1 

n  24  yV2  .  Q 

SrQ  =  - — -  c  sin0  . 

This  solution  is  used  to  check  the  accuracy  of  the  numeri¬ 
cal  solution. 

3 . 2  Axi-Symme tr ic  Flow  in  a  Cone 

The  second  problem  considered  is  the  axi-symmetr ic , 
quasi-static  flow  of  an  incompressible  material  in  a  cone. 
Again,  the  flow  is  assumed  to  be  radial.  Spherical  coordinates 
r,  <J>  ,  and  0  (Figure  3.2),  are  used. 

The  non  zero  components  of  the  rate  of  strain  tensor, 
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FLOW  IN  A  CONE 


in  spherical  coordinates,  are 


d 


r 


9vr 

3r~ 


d<j>  =  de  =  ■—  »  0.12) 

and 

tcf)  2r  3  <t> 


From  the  continuity  equation, 


' 
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3  v 


r 


+  2 


the  radial  velocity  is  found  to  be  of  the  form 


vr 


(3.13) 


The  components  of  the  rate  of  strain  tensor,  equations  (3.12), 
become 


d 


r 


and 


d0 


 g  (<f>  ) 


g  ’  (4>) 

2r3  ’ 


(3.14) 


upon  substitution  of  relation  (3.13).  The  prime  denotes 
differentiation  with  respect  to  <f)  . 

The  equations  of  motion  are 


3ar  i  9  ar<p  +  1^ 

3  r  r  34>  r 


(2  ar 


a0  +  arcf>  cot  ^ 


0  , 


and 


(3.15) 


which  become,  with  the  substitution  (2.2)  and  noting  that 


;0  =  -  (sr  +  s(f>) 


3r  3  r 


3  S  S 

+  3  ~  + 


3  p  _  o  a  r  ,  0  °r  ,  1 


3  S  r(f)  1 

F  ■  r  -?T  +  r  COt  *  > 


and 


IP  =  r  i!li  +  i!l  +  3S 

3*  3  r  3  * 


r  (p  ' 


(3.16) 


Elimination  of  the  hydrostatic  pressure  yields 


a2sr  +  2  dsr  _  d2s<p  +  i 
3r34>  r  9<J>  3r3(J)  r  9^2 


+  i 


-  r 


^ _ ^rcj) 

3r2 


+  1  ±  <Sr*  cot*)  -4^=0, 


(3.17) 


which,  again,  is  entirely  in  terms  of  the  deviatoric  stresses 
Substituting  equations  (3.14)  into  the  power  law  gives 


5-n 

S  =  - 
S  r  „  3n 


n-1 


yg [ 12g2  +  (g 1 ) 2  ]  2  , 


3  -n 


n-1 


=  S0  =  Pg[12g2  +  (g')2]  2  , 
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and 


*  r  4> 


1-n 
2  2 
^  3n 


n-1 


Pg  '  [ 12g2  +  (g ' ) 2 ]  2 


where  Sr  ,  ,  Sq  ,  and  Sj.^  are  the  non-zero  components  of 

the  deviatoric  stress  tensor,  and  the  prime  denotes  differen¬ 
tiation  with  respect  to  <j>  . 

Letting 

m  (  <J> )  =  8  (  ft) .  , 

V 

where  m  is  a  non-dimensional  velocity,  and  V  is  a  constant 
to  be  determined  from  the  volume  flow  rate,  the  expressions 
for  the  stress  deviators  become 

5-n 

sr  ■  73 r- yV"  mG- 

3-n 
?  2 

S(f)  =  Se  =  ~~~  yVn  mG  ,  (3.18) 


and 


1-n 


>r  (p 


yVn  m'G 


3n 


9 


where 


n-1 

G  =  [12m2  +  (m')2]  2 
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Substitution  of  equations  (3.18)  into  equation  (3.17) 


gives 


(m'G)"  +  (ra'G  cot  <J> )  '  +  9n(l-n)  m  ’  G 

+  6(3n-2)  (mG) '  =  0  ,  (3.19) 

and  the  expression  for  the  hydrostatic  pressure  becomes 

1-n 

p  =  2  2.  uVn  [—  (1-n)  mG  -  {(m'G)'  +  m'G  cot  4> }  ]  +  A.  (3.20) 

F  r3n  K  n  3n 

Equations  (3.19)  and  (3.20)  reduce  to  those  for  the 
quasi-static  flow  of  a  Newtonian  fluid  in  a  cone  for  n  -  1, 
and  De  Vries  (18)  showed  that  they  reduce  to  the  equations 
obtained  by  Shield  (8)  for  the  r ig id-p er f ec t ly  von  Mises 
solid  for  n  =  0.  Dividing  equation  (3.18)  by  3n,  integrat¬ 
ing,  and  rearranging,  yields 


-  (1-n)  mG  -  {(m'G)’  +  m'G  cot  <|>  }  m 

n  3n 

=  2mG  +  3 ( 1 -n )  /  m'G  d$  +  A. 

0 


Substituting  this  expression  in  equation  (3.20)  and  evaluat 
ing  at  n  =  0  gives 


2  7 


P 


/2  y[2mG  +  3  /  m’G  dc|>  ]  +  E(r)  +  A  . 

0 


This  is  equivalent  to  the  expression  for  the  pressure 
from  Shield’s  equations,  since,  for  n  =  0, 

mG  =  +  -i—  /l-(m'G)2  , 

“  2/3 

and  with  the  substitution 

T  =  m’G  , 

the  expression  for  the  pressure  becomes 


_  4) 

p  =  /2  y  [  -  — L  /l-x2  +  3  /  T  d<|)  ]  +  E(r)  +  A  . 

/  3  o 


The  function  E(r)  can  be  evaluated  since,  from 
tion  (3.20), 


(l£)  =  [  (m '  G )  ’  +  m’G  cot  <j>  -  12mG], 

9r  n=  0  r 


and  from  equation  (3.19) 


(m*G)  '  +  m'G  cot  4>  -  l2mG  =  0  . 


f  ound 


(3.21) 


equa- 


From  equation  (3.21) 
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(f£)  =  E’(r)  , 

3r  n=0 


which  implies  that 


E  ( r )  =  9/2  y  D  In  r  +  B  , 


and  the  expression  for  the  pressure  becomes 


-  /2  „[- 


T  2  _|_ 


/ 

0 


T  d <f>  +  c  In  r  ]  +  A 


The  problem  is  reduced  to  solving  equation  (3.19) 
subject  to  the  boundary  conditions 

m  ( 0 )  =  1  , 

m'  (0)  =  0 

from  symmetry,  and  the  no  slip  boundary  condition 

m (a)  =  0  , 

where  a  is  the  semi-angle  of  the  cone.  Again,  in  general, 
this  must  be  done  numerically. 

2 

De  Vries  (18)  showed  that,  for  n  =  -j ,  a  closed  form 
solution  of  equation  (3.19)  in  terms  of  m*  G  can  be  obtained. 


The  solution  is 
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m '  G  =  B  sincf)  , 


and  the  shear  stress  becomes 


S 


r  <J) 


1  2 
26  yV3 

r3 


B 


s  in<j> 


which  is  used  to  check  the  accuracy  of  the  numerical  solution. 


CHAPTER  IV 


NUMERICAL  METHOD 

In  the  preceding  chapter,  the  governing  equations  were 
derived  for  the  flow  problems  under  consideration.  These 
equations  can  be  integrated  numerically  using  Gill’s 
variation  of  the  Runge-Kutta  fourth-order  method.  The 
problem  must  first  be  converted  to  an  initial  value  problem 
by  assuming  a  value  for  the  unknown  initial  condition  m"(o) . 

A  trial  and  error  procedure  to  find  the  value  of  m"(o)  which 
satisfies  the  no-slip  boundary  condition  is  very  costly  in 
computer  time.  A  method  developed  by  Clutter  and  Smith  (19) 
was  adapted  to  the  present  problem  to  systematize  the  solution 
of  these  equations. 

4 . 1  Reformulation  as  an  Initial  Value  Problem 

The  Runge-Kutta  method  can  be  used  to  solve  a  system 
of  first  order  ordinary  differential  equations  with  pre¬ 
scribed  initial  conditions.  The  two  problems  under  consider¬ 
ation  must  first  be  put  in  this  form. 

The  plane  flow  problem,  equation  (3.9)  and  the  asso¬ 
ciated  boundary  conditions,  yields  the  system  of  first  order 
equations 


y  i '  (9)  =  y2  )  > 


y  2  '  ( 9 )  =  y  3  (9  )  > 
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and 


y  3  T  ( 0  )  =  F  [y-j^  (0  )  ,y  2  (6  )  ,y  3  (0  )  ]  > 


where 


y  ^  (0  )  =  m(0)  , 


y  2  ( 0  )  =  tn  *  (0  )  , 


y  3  (  0 )  =  m"(0)  , 


F[yi(9)  ,y2(e)  ,y3(9)]  -  - - - \  {(4n2 -12n+4)  m'M 

M+ (n-1) (m  ) 


-  [(4n^-6n+2)  m  +  (n-l)m"]  M' 


-  (n-1) (n-3)  m' (M' ) 2 / 4M 


-  (n-1)  m’[4(m')2  +  4mm"  +  (m")2]}  , 

and 

M  =  4m^  +  (m ' ) ^  . 


The  initial  values  are 


y  1  ( 0 )  =  1  , 
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y 2 (°)  =  0  > 


and 


y 3 (0)  =  3  • 

3  is  the  unknown  value  of  m’(0)  to  be  determined  from  the 
no  slip  boundary  condition. 

Similarly,  for  the  axi-symmetr ic  flow  problem,  equation 
(3.18)  and  the  associated  boundary  conditions,  gives  the 
system  of  first  order  equations 

Yl  ’  ($)  =  Y2  > 

Y2  '  (4>)  =  Y 3  (4>  )  » 


and 


y 3 1  (4> )  =  F[<f),yi(<t>) » y 2 ( )  » y 3 ^ )  ]  » 


wher  e 


f [<t> , yi (<i>)  ,y 2 > Y2 i  ~  — 7 — 77 — mr  ^9n2  27n+12^  m' 

M+ (n-l)(m  ) 


+ 


CSC 


(J)  m '  ]  M 


cot  ({)  m" 
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-  [(9n2-15n+6)  m  +  (n-1)  m" 


+  cot  4>  m'  ]  M 1 


-  (n-1) (n-3)  m’(M')^/4M 

-  (n-1)  m,[l2(m’)2  +  12mm"  +  (m")2]), 
M  =  12m2  +  (m' )2 , 


y1(<|))  =  m(<J>)  , 

y2(4>)  =  m’(<j>)  , 


and 


y  3  ( )  =  m"((j))  , 


The  initial  conditions  are 


y  1  ( 4> )  =  l  > 


y2(4>)  =  0  , 


and 
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y3(<J>)  =  3  . 

Again,  3  is  to  be  determined  from  the  no  slip  boundary 
condition . 

4 

4 . 2  Determination  of  the  Unknown  Initial  Condition 

A  purely  trial  and  error  procedure  to  satisfy  the  no 
slip  boundary  condition  uses  too  much  computer  time.  Clutter 
and  Smith  (19)  encountered  the  same  difficulty  in  solving 
the  compressible  boundary  -  layer  equations  for  an  arbitrary 
axi- s y mme t r i c  body.  The  problems  considered  here  are 
simpler  since  the  flow  is  incompressible  and  thermal  effects 
are  neglected.  The  method  can  be  applied  to  both  the  plane 
and  axi-symme tr ic  flow  problems. 

The  procedure  is  as  follows.  A  value  for  3  is  first 
assumed,  and  the  system  of  equations  integrated.  The  initial 
value  for  3  is  increased  if  m(a)  is  greater  than  zero,  or 
decreased  if  m(a)  is  less  than  zero.  The  equations  are  then 
integrated  again  using  the  new  value  for  3.  When  both  a 
high  and  a  low  value  for  m(a)  is  obtained,  the  two  values  for 
3  are  averaged,  and  the  equations  integrated  again  using  this 
averaged  value.  This  averaging  process  could  be  continued 
until  the  magnitude  of  m(a)  is  as  small  as  desired.  This 
process  is  speeded  up  considerably  by  using  a  three  point 
interpolation  formula.  When  three  solutions,  m^(0), 
and  m^ ( 6 ) ,  are  obtained  such  that 

-  k_<m^(a)  _<  k  , 


35 


-  k  <_  (ot)  ^  k  , 


and 


-  k  £  m3(a)  <_  k  , 

where  k  is  a  prescribed  limit,  and  two  of  m^(a),  m2  (oO  ,  and 
m3  (a)  are  greater  than  zero  and  the  third  less  than  zero,  or 
conversly,  they  are  combined  by  an  interpolation  formula 
which  automatically  satisfies  all  the  boundary  conditions. 
The  interpolated  solution  can  be  checked  by  integrating  the 
equations  using  the  interpolated  value  for  3.  Clutter  and 
Smith  found  the  computer  time  required  using  this  method  to 
be  about  one  fourth  that  required  if  the  averaging  process 
wa s  used. 

The  interpolated  solution  is  found  by  substituting  the 
three  solutions  in  the  formulas 

m(0)  =  A^m-^(0)  +  A2fn2(9)  +  , 

m  ’  (0  )  =  A3_m]_'(0)  +  A2m2  '  (9  )  +  A3m3*(0)  , 


and 


m" (0 )  =  A1m1"(0)  +  A2m2"(e)  +  A3m3"(0)  , 


wher  e 
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and 


Ai  ■ 


m2(ot)  1113(a) 


A  9  - 


An  = 


[m^a)  -  m2(a)  ]  [m^(a) 

-  m3 (a) ] 

mi(a)  1113(a) 

[ tn 2  (a)  -  (a)  ]  [m2  (a) 

-  m3 (a)  ] 

m^  (a )  m2  (a ) 

3  [m3  (a)  -  (a)  ]  [m3  (a)  -  m2  (a.)  ]  * 


By  substitution,  it  is  easy  to  see  that  all  the  boundary 
conditions  are  satisfied  by  the  interpolated  solution. 

The  other  quantities  of  interest,  pressure  and  shear 
stress,  can  be  found  in  the  same  way. 

The  value  of  k  must  be  chosen  to  give  the  required 
accuracy.  A  value  for  k  of  0.2  was  found  to  give  solutions 
which  agreed  to  four  decimal  places  with  the  solution 
found  by  integrating  the  equations  using  the  interpolated 
value  for  3.  Also,  good  agreement  with  the  known  solutions 
discussed  in  the  preceding  chapters  was  obtained.  Some 
experimentation  was  necessary  to  determine  the  value  of  k 
which  gave  sufficiently  accurate  solutions  using  a  minimum 
of  computer  time. 

The  value  of  3  assumed  initially  is  important  in 


reducing  computer  time.  Since  the  solutions  for  both  plane 
flow  in  a  converging  channel  and  a xi-symme tr ic  flow  in  a  cone 
are  known  for  a  Newtonian  fluid,  for  which  n  =  1,  and  it  was 


' 
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desired  to  solve  the  equations  for  a  series  of  values  for  n 
less  than  one,  the  known  value  for  3  for  n  =  1  was  used  as 
the  first  trial  value  for  3  for  n  =  0.9.  A  series  of 
solutions  for  n  =  1  down  to  n  =  0.01,  in  steps  of  0.1,  was 
obtained  in  this  way.  Only  four  trials  were  necessary  to 
obtain  the  three  solutions  needed  for  use  in  the  inter¬ 
polation  formula  for  each  value  of  n. 

Some  trial  solutions  and  a  comparison  of  the  inter¬ 
polated  solution  to  the  solution  found  by  integrating  the 
equations  using  the  interpolated  value  for  3,  for  plane  flow, 
a  semi-angle  of  fifteen  degrees,  and  n  =  0.5,  are  presented 
in  the  next  chapter. 

4 . 3  The  Special  Case:  n  =  0 

No  difficulty  is  encountered  in  obtaining  solutions 
for  arbitrarily  small  values  of  n.  However,  the  computer 
time  needed  increases  with  decreasing  n,  since  smaller 
step  sizes  are  required.  When  n  =  0  the  situation  is  differ¬ 
ent.  The  velocity  derivative  is  infinite  at  the  boundary, 
and  therefore  the  numerical  solution  can  not  be  extended 
right  to  the  boundary.  Two  approaches  are  possible. 

1.  The  zero  velocity  boundary  condition  can  be 
satisfied  at  a  point  arbitrarily  close  to  the 
boundary . 

2.  The  values  of  3  satisfying  the  zero  velocity 
boundary  condition  for  n  £  1  can  be  extrapolated 
to  zero,  and  the  extrapolated  value  used  to 
integrate  the  equations  for  n  =  0. 

The  first  approach  is  unattractive  since  too  much  computer 

time  is  required  to  get  a  solution  sufficiently  close  to  the 


, 
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boundary.  The  second  approach  was  taken  by  De  Vries  (20). 

A  shear  stress  boundary  condition  can  also  be  imposed 
and  the  same  numerical  methods  discussed  above  employed. 

Good  agreement  with  Nadai's  solution  for  the  r i g i d -per f ec t ly 
plastic  Von  Mises  solid  was  obtained  in  this  way  for  a  semi¬ 
angle  of  15  degrees. 


CHAPTER  V 


RESULTS 


The  results  of  the  numerical  analyses  are  presented 
in  this  chapter.  First,  some  tables  are  given  to  demonstrate 
the  merit  of  the  numerical  method,  and  second,  the  solutions 
obtained  are  given  in  the  form  of  graphs  for  semi-angles 
of  5,  10,  and  15  degrees,  and  various  values  of  n. 

Table  5.1  gives  the  trial  values  for  3  used  to  obtain 
the  solution  to  the  plane  flow  problem  for  a  semi-angle  of 
15  degrees  and  n  equal  to  0.5.  The  first  trial  value  is 
the  value  of  3  which  satisfies  the  no  slip  boundary  condition 
for  n  equal  to  0.6.  The  next  three  trials  gave  the  solutions 
required  for  substitution  in  the  interpolation  formula.  As 
a  check  on  the  validity  of  the  interpolated  solution,  the 
equation  was  integrated  using  the  interpolated  value  for  3. 
The  solutions  obtained  by  the  two  methods  are  compared  in 
table  5.2.  The  solutions  agree  to  four  significant  figures. 

A  further  check  on  the  method  is  to  compare  the  numerical 
solutions  obtained  for  the  ax i-symme tr ic  flow  problem  for 
n  equal  to  JL,  and  for  the  plane  flow  problem  for  n  equal  to 


—  ,  to  the  analytic  solutions.  This  is  done  in  table  5.3, 
and  again  there  is  satisfactory  agreement.  These  results 
show  the  validity  of  the  method.  Closer  agreement  could  be 
obtained  by  decreasing  the  value  of  k,  with  a  consequent 


increase  in  computer  time. 

Figures  5.1  through  5.9  present  the  results  obtained 
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for  plane  flow  in  a  converging  channel  for  semi-angles  of 
5,  10  and  15°,  and  figures  5.10  through  5.15  the  results 
for  axi-symmetr ic  flow  in  a  cone,  for  semi— angles  of  5  and 
10°,  and  various  values  of  n. 

Figures  5.16  through  5.19  give  the  results  for  the 
plane  flow  of  a  r igid-per f ectly  plastic  Von  Mises  solid  in  a 
converging  channel  with  a  shear  stress  boundary  condition. 
The  results  agree  with  Nadai ' s  solution. 
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TABLE  5.1 

TRIAL  SOLUTIONS  FOR  PLANE  FLOW 
n  =  0.5  a  =  15° 


Trial 

1  -19 

2  -16 

3  -18 

4  -17 

Interpolated  -17 


3 

v r  (a) 

.  9814 

-0 . 245150 

.9814 

+0.073443 

.4814 

-0 . 078910 

.6909 

+  0 .003114 

.7211 

0 . 000000 

-  ' 
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TABLE  5.2 

COMPARISON  OF  INTERPOLATED  AND 
INTEGRATED  SOLUTIONS  FOR  PLANE  FLOW 

a  =  15° 


0  m  (0  ) 


Degrees 

Interpolated 

Solution 

Integrated 

Solution 

1 

1.000000 

1.000000 

3 

0.975194 

0.975193 

6 

0.893993 

0.893990 

9 

0.732847 

0 .732843 

12 

0.450761 

0.450757 

15 

0.000000 

P  (6) 

0.000003 

1 

0 . 000000 

0 . 000000 

3 

0.130321 

0 . 130323 

6 

0.501925 

0 . 501933 

9 

1.028350 

1.028397 

12 

1.599365 

1.599467 

15 

2 . 189448 

S  n 
r  0 

2.189571 

1 

0.000000 

0 . 000000 

3 

0.779888 

0 .779913 

6 

1.557638 

1 .557689 

9 

2.331119 

2.331195 

12 

3 . 098210 

3.098312 

15 

3.856809 

3.856936 

4  3 


TABLE  5.3 

COMPARISON  OF  THEORETICAL  AND 
NUMERICAL  SOLUTIONS 

a  =  15° 


Axi- 

Symmetric  Flow 

(n=2 / 3 ) 

S  A 
r  (p 

Degrees 

Theoretical 

Numerical 

0 

0 

1 

0 . 32179 

0 .32180 

2 

0 . 64349 

0 . 64350 

3 

0 . 96498 

0 .96500 

4 

1 .28619 

1 . 28621 

5 

1 . 60700 

1 .60703 

6 

1 .92732 

1 .92736 

7 

2 . 24706 

2 .24711 

8 

2 .56611 

2 .56616 

9 

2.88438 

2.88444 

10 

3 .20177 

3 . 20184 

11 

3 .51818 

3 .51826 

12 

3.83353 

3 . 83361 

13 

4.14770 

4.14779 

14 

4.46062 

4 .46071 

15 

4.77217 

4 .77227 
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TABLE  5.3  (CONTINUED) 

COMPARISON  OF  THEORETICAL 

AND 

e 

NUMERICAL  SOLUTIONS 

a  =  15 

Plane  Flow  (n=l/2) 

S  A 

r  0 

;rees 

Theoretical 

Numerical 

0 

0 

0 

1 

0 . 26008 

0 .26007 

2 

0.52007 

0 .52006 

3 

0.77991 

0.77989 

4 

1.03952 

1.03948 

5 

1.29880 

1.29876 

6 

1.55769 

1.55764 

7 

1.81610 

1.81604 

8 

2.07397 

2.07390 

9 

2.33120 

2.33112 

10 

2.58772 

2 . 58763 

11 

2.84345 

2.84335 

12 

3.09831 

3 . 09821 

13 

3.35223 

3.35212 

14 

3.60513 

3.60502 

15 

3.85694 

3.85681 

RESULTS 


FOR 


PLANE  FLOW 


m{0) 
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FIGURE  5.1 

tn  (0 )  VERSUS  0  :  a  =  15* 


m(0) 
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FIGURE  5.2 

m  (0)  VERSUS  0  :  a  =  10° 


m(e) 
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FIGURE  5.3 

m  (0)  VERSUS  0:  a  =  5° 
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FIGURE  5.4 
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-  s  A (0)  VERSUS  0 :  a  =  15 
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FIGURE  5.5 
2n 

- -  S  _  ( 0 )  VERSUS  0:  a  =  10 

„n  r  0 
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n  =  0.04 
n  =  0.02 
n  s  0.09 
n  =  0.0  1 


0 


FIGURE  5.6 


2n 


S  Q  (0 )  VERSUS  0:  a  =  5‘ 


,n  r  0 
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FIGURE  5.7 
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FIGURE  5.8 
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FIGURE  5.9 
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AXI- SYMMETRIC  FLOW 
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FIGURE  5.10 

m  (0 )  VERSUS  0 :  a  =  10" 
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FIGURE  5.11 


5° 


m  (  0 )  VERSUS  0: 


a 


Sre(e) 
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e 

FIGURE  5.12 
3  n 

- -  S  Q(0)  VERSUS  0:  a  =  10" 

n  r0  v 

yv 


Sr  9(e) 
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FIGURE  5.13 
3n 

I -  S  „  (0 )  VERSUS  0:  a  =  5° 

yVn  r0 
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PLANE  FLOW 
n  =  0 
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FIGURE  5.16 

SHEAR  STRESS  BOUNDARY  CONDITION 
3  VERSUS  n :  a  =  15  ,  n  =  0 


64 


FIGURE  5.17 

SHEAR  STRESS  BOUNDARY  CONDITION 
m(0)  VERSUS  0  :  a  =  15'  ,  n  =  0 


65 
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FIGURE  5.18 

SHEAR  STRESS  BOUNDARY  CONDITION 


S  a  t0) 

r  6 
U 


VERSU  S  0 :  a  =  15 
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Q> 

a 


e 


FIGURE  5.19 

SHEAR  STRESS  BOUNDARY  CONDITION 
P VERSUS  6:  a  =  15',  n  =  0 
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CHAPTER  VI 


INLET  EFFECTS 


In  order  to  obtain  solutions  to  the  pseudo-plastic 
flow  problems  considered  in  this  thesis,  it  was  necessary 
to  assume,  from  the  outset,  fully  developed  flow.  That 
the  streamlines  for  fully  developed  flow  are  radial  was 
verified,  since  solutions  were  obtained.  The  velocity 
profiles  given  by  these  solutions  are  established  at  some 
distance  from  the  entrance,  but  to  find  this  distance  it 
would  be  necessary  to  consider  the  complete  equations  of 
motion  and  the  continuity  equation.  This  is  a  system  of 
three  coupled,  second  order,  non-linear  partial  differential 
equations  which  are  very  complex.  An  attempt  was  made  to 
solve  a  simplified  form  of  these  equations  for  axi-symmetr ic 
flow  of  a  Newtonian  fluid  in  a  cone  by  a  method  suggested  by 
Sutterby  (21).  This  method  appears  to  be  inadequate. 

Another  method,  based  on  variational  principles,  offers  some 
hope  for  obtaining  approximate  solutions  for  entrance  effects. 
This  is  an  area  for  future  work  in  this  field. 

6 . 1  The  Newtonian  Fluid 

Both  Sutterby  (21)  and  Tanner  (13)  have  investigated 
entrance  effects  for  the  flow  of  a  Newtonian  fluid  in  a 
cone.  Sutterby  applied  a  numerical  technique,  while  Tanner 
used  analytic  methods.  In  both  cases  the  results  obtained 
appear  to  be  neither  conclusive  nor  very  satisfactory. 
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For  small  semi-angles,  Sutterby  assumed  that  the 
tangential  velocity  would  be  small,  and  reduced  the  equations 
of  motion  and  continuity  to 


Pv, 


3  Vi 

3  r 


-  ±E  + 

dr 


y  rl  3  /n  1 

,2  0  39  (6  36  )] 


1  _3_ 
r  3  r 


+  9  30  <0V0) 


0  , 


and 


a 

Q  =  2irr^  J  vr  0  d0  . 
0 


For  entrance  conditions  a  constant  pressure,  a  zero  tangen¬ 
tial  velocity,  and  a  parabolic  profile  for  the  radial  velo¬ 
city, 


vr  =  V[1  -  (£)2]  , 

were  assumed.  For  small  semi-angles  this  is  very  close  to 
the  radial  flow  solution  given  by  equation  (1.6).  Using 
finite  differences,  the  equations  were  solved  for  various 
values  of  the  Reynolds  number.  The  results  showed  that,  as 
the  Reynolds  number  approaches  zero,  the  solution  approaches 
the  quasi-static  radial  flow  solution,  and  for  large 
Reynolds  numbers  it  approaches  the  ideal  inviscid  radial 
flow  solution.  However,  a  solution  for  a  Reynolds  number  of 


' 
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zero  was  not  given;  it  was  assumed  that  the  known  radial  flow 
solution  would  be  the  solution,  regardless  of  conditions  at 
the  entrance  and  the  neglect  of  the  variation  in  the  pressure 
with  6.  This  method  was  found  to  yield  results  for  the 
parabolic  inlet  velocity  profile,  and  for  an  arbitrarily 
small  Reynolds  number,  which  appeared  to  be  correct.  However, 
for  zero  Reynolds  number,  the  solution  did  not  approach  the 
radial  flow  solution  at  some  distance  from  the  entrance  as 
was  expected.  Inclusion  of  the  variation  of  pressure  in 
the  0  direction  had  an  insignificant  effect  on  the  results. 

A  more  severe  test  of  the  method  was  tried  using  the  inlet 
velocity  profile 

vr  =  V [ 1  -  <£)10]  . 

This  change  in  the  inlet  conditions  was  not  reflected  in  the 
results.  This  led  to  the  conclusion  that  the  equations  used 
by  Sutterby  are  over-simplified,  and  are  inadequate  to 
describe  the  flow,  at  least  in  the  low  Reynolds  number 
regime,  and  that  reasonable  solutions  were  obtained  only 
because  of  the  choice  of  a  parabolic  inlet  velocity  profile. 

To  find  the  distance  from  the  inlet  at  which  the  flow  is  fully 
developed,  a  much  more  complicated  form  of  the  equations 
must  be  solved  . 

Tanner  (13)  attempted  to  find  the  effect  of  entrance 
conditions  by  introducing  Stokes'  stream  function  in  the 
inertialess  equations  of  motion  in  spherical  coordinates. 


which  yields 
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s in0 

9  r 


( 


s  inO 


sin0  9  ✓  1  9  x  ,  2 

r 2  90  sin0  90; J 


* 


0  , 


where 


v  =  .  i  2! 

r  r^sine  30  ’ 

and 


vp--  l—_  21 

®  r  si n(9  9r 


The  solution 


^  =  I  rm  fm(x) » 
m 


where 

X  =  cos0  , 

was  assumed,  and  solutions  for  the  fm  in  terms  of  Legendre 
functions  were  obtained.  The  values  of  m  which  satisfy  the 
boundary  conditions  were  found  to  be  complex.  These  values 
of  m  were  used  to  find  the  distance  at  which  inlet  effects 
would  no  longer  be  felt.  For  semi-angles  smaller  than 
thirty  degrees.  Tanner  found  this  distance  to  be  approximately 
equal  to  the  radius  of  the  cone  at  the  inlet,  but  the  manner 
in  which  this  conclusion  was  reached  is  not  clear. 

The  problem  of  entrance  effects  for  the  quasi-static 


?1B5  1  -  -  - ,  , . 
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flow  of  a  Newtonian  fluid  in  a  cone  cannot  be  considered 
solved..  Experimental  work  and  refined  numerical  analysis 
are  needed  to  settle  the  question. 

6 . 2  Variational  Principles 

The  problem  of  entrance  effects  for  a  Newtonian  fluid 
is  complex,  and  satisfactory  solutions  have  not  yet  been 
obtained  for  the  flows  considered  in  this  thesis.  The 
equations  describing  the  flow  of  non-Newtonian  materials, 
and,  in  particular,  power  law  materials,  are  much  more 
complicated,  and  there  seems  to  be  little  hope  for  obtaining 
exact  solutions  for  entrance  effects. 

Bird  (22)  showed  that  the  equations  of  motion  and 
continuity  for  the  quasi-static  flow  of  an  incompressible 
power  law  material  are  equivalent  to  the  Euler -Lagr ang e 
equations  for  the  integral 


Iff  F  dv  , 


(6.1) 


where 


F 


*2 

(  n  + 1 ) 


(p  +  p(f>)  I  , 


n-1 

n  =  2y  (2l2 )  2 


4>  is  the  potential  function  for  the  body  forces,  and  I  ^ 


and 


. 
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I  are  the  first  and  second  invariants  of  the  rate  of  strain 
tensor.  The  velocity  must  be  specified  over  the  entire 
surface  of  the  volume  over  which  (6.1)  is  taken. 

With  this  variational  principle  approximate  solutions 
may  be  obtainable.  If  a  velocity  distribution,  containing 
arbitrary  constants,  that  satisfies  the  no-slip  boundary 
condition  and  gives  the  prescribed  velocity  profile  at  the 
inlet,  and  fully  developed  flow  at  some  distance  down  stream, 
can  be  constructed,  then  the  arbitrary  constants  may  be 
chosen  to  minimize  (6.1).  Even  for  simple  velocity  dis¬ 
tributions  the  expression  for  F  will  be  very  complicated, 
and  even  if  a  solution  can  be  found  in  this  way,  there  is 
no  way  of  knowing  how  close  the  approximate  solution  is  to 
the  exact  solution.  Convergence  could  be  checked  by  taking 
more  arbitrary  constants  in  the  assumed  velocity  distribution. 
This  variational  principle  also  opens  the  possibility  of 
using  the  finite  element  technique.  This  is  an  area  for 
further  investigation. 


CHAPTER  VII 


CONCLUSION 


In  this  thesis  the  power  law  was  considered  as  an 
empirical  approximation  for  the  constitutive  equations  of 
pseudo-plastic  materials  for  a  certain  range  of  flow  rate. 
Experimental  results  are  needed  to  compare  with  the 
theoretical  results  obtained  in  order  to  check  the  validity 
of  the  power  law.  The  assumption  that  the  streamlines  are 
radial  for  fully  developed  flow  was  justified  since  solutions 
were  obtained,  but  the  problem  of  entrance  and  exit  effects 
still  remains  unsolved.  More  work  needs  to  be  done  in  this 
area  . 

The  numerical  technique  developed  in  this  thesis  can  be 
adapted  to  most  two  point  boundary  value  problems  for  ordinary 
differential  equations.  It  appears  to  be  more  straight  forward 
and  more  universally app li cab le  than  other  techniques,  such  as 
quasi-linear iza tion  and  invariant  imbedding,  and  is  a  great 
deal  more  economical  in  computer  time  than  the  usual  trial 
and  error  procedure.  Hence,  it  should  prove  of  value  in 
solving  other  two  point  boundary  value  problems. 


■ 
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APPENDIX 


APPENDIX 


AN  EXAMPLE  COMPUTER  PROGRAM 

The  program  developed  to  obtain  the  numerical  solution 
to  the  plane  flow  problem  for  a  semi-angle  of  15  degrees  is 
presented  below  for  reference  purposes. 


o  o  o  o  o  o 
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DOUBLE  PRECISION  PRMT(5),  Y(3),  D ERY ( 3 )  ,  AUX  (8,3), 
BETA,  ALPHA, 

1TN,  Yl,  YY1 ,  YYY1 ,  YYYYl,  BETL ,  BETH,  DEL,  Vl(31), 

V2  (31) ,  V3 ( 31) 

2,  PI (31) ,  T1  ( 31)  ,  Al,  A2 ,  A3,  PO,  Y2(3,31),  YY(3,31), 

YYY(3 ,31)  , 

3P ( 3 , 31)  ,  T ( 3 , 31) 

EXTERNAL  FCT ,  OUTP 

COMMON  TN  ,  BETA  ,  Y 1 ,  PO , Y 2 , YY , YYY , P , T , LL , LLL 
ND IM=  3 


SET  PARAMETERS 
TEST=0 . 2DO 


ND  =  15 
NL=  2*ND+1 


NN=ND* 10 
ALPHA=0 . 15D2 

CALCULATE  BETA  FOR  A  NEWTONIAN  FLUID 

BETA=0 .  4D  1/ (DCOS (2 .* ALPHA* 3 .14159  265  35  89  79  3/180  .) -0  .ID  1) 


C 


C 


C 

C 

C 


PREPARE  FOR  DRKGS 

SET  VALUES  OF  PRMT  AND  DERY 

PRMT (1) =0 . 0D0 

PRMT (2) = ALPHA* 3 .141592653589793/180. 
PRMT (3) =PRMT (2) /FLOAT (NN) 

PRMT (4)=0 . ID -5 

SET  VALUE  OF  EXPONENT 

DO  28  JK=10 , 100 , 90 

DO  28  IN=1 , 10 

TN=FLOAT (11-IN) /FLOAT ( JK) 

JC  1  =  1 
JC  2=1 
LLL=  1 


DEL=0 . 3D1 
BETL=0 . 0D0 
BETH=0 . 0D0 
WRITE(6 ,100) 

WRITE (6 ,101)  ALPHA, TN 
WRITE ( 6 ,108) 

10  CONTINUE 
LL  =  0 

SET  INITIAL  VALUES 


Y(1)=0.1D1 

Y(2)=0.0D0 

Y  (  3 —  BETA 

PO=2**((l.-TN)/2.)*(2./TN*(l.-TN)*Y(l)-l./(2.*TN)* 

1(Y(3)+Y(2)*(TN-1.)/2./(4/*Y(1)**2+Y(2)**2)*(8.*Y(1) 

2*Y (2) +2 .*Y(2)*Y(3) ) ) ) * (4 .*Y (1) **2+Y (2) **2) ** ( (TN-1 . 

DO  11  J=l, 3 
11  D  ERY (J)=0.1D1/0.3D1 


/  2  .  ) 


CALL  DRKGS (PRMT , Y ,DERY , DNIM , IHLF ,FCT ,OUTP , AUX) 


TERMINATE  PROCEDURE  IF  BOUNDARY  CONDITION  IS  SATISFIED 


o  n 
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C 

C 


IF  (DABS (Yl)  .LT . 0  .ID -4)  GO  TO  27 


BEGIN  PROCEDURE  TO  FIND  VALUES  FOR  USE  IN  INTERPOLA¬ 
TION  FORMULA  ,  x  .  nrv  os 

IF (DABS  (BETL)  .GT .0 .ID -8 .AND .DABS (BETH)  . GT .0 .ID -8) 


GO  TO  15 

IF(DABS(Y1) .LT.TEST)  GO  TO  16 

12  IF  (Yl )  13,27,14 

13  BETL=BETA  , 

IF (DABS (BETL) . GT . 0 . ID -8 .AND .DABS (BETH) 

GO  TO  15 
BETA=BETA+DEL 


.GT.0.1D-8) 


GO  TO  10 


14 


15 


BETH=BETA 

IF  (DABS (BETL)  . GT .0 .ID -8 .AND .DABS 
GO  TO  15 
BETA=BETA-D  EL 
GO  TO  10 

BETA= (BETL+BETH) /0 . 2D  1 
BETL  =  0 . 0D0 
BETH=0 . OD  0 
DEL=DEL/0 . 3795D1 


(BETH) .GT . 0 . ID -8 ) 


GO  TO  10 

16  IF(Y1)17,27 ,22 

17  IF (JC1+JC2 ,EQ  .5)  GO  TO  25 
IF ( JC1 . LT . 2)  GO  TO  18 

IF ( JCl . LT . 3)  GO  TO  20 
GO  TO  12 

18  YY1=Y 1 

IF (LLL  .  LT3)  LLL=LLL+1 
JC 1= JC 1+1 

IF (JC1+JC2  .EQ  .5)  GO  TO  25 

20  IF (DABS (Yl -YY 1 )  . LT . 0 . ID -11 . OR . DABS ( Yl-YYYY 1 )  .LT.0.1D-11) 

GO  TO  12 
YYY1=Y1 

IF  (LLL. LT. 3)  LLL=LLL+1 
IF (Yl .LT . 0 . 00D0)  JCl=JCl+l 
IF  (Yl .  GT  .0 .00D0)  JC  2  =  JC  2+1 
IF  (JC1+JC2 .EQ . 5)  GO  TO  25 
GO  TO  12 

22  IF (JC1+JC2 .EQ .5)  GO  TO  25 
IF  ( JC 2  . LT . 2 )  GO  TO  23 

IF  ( JC2  . LT . 3)  GO  TO  20 
GO  TO  12 

23  YYYY1=Y 1 

IF (LLL. LT. 3)  LLL=LLL+1 

IF ( JC1+JC 2 . EQ . 5 )  GO  TO  25 
GO  TO  12 


25 


CALCULATE  COEFFICIENTS 
A1=Y2(2,NL)*Y2(3,NL) / ( 
Y2  (3 ,NL) ) ) 


FOR  INTERPOLATION  FORMULA 
(Y2 (1,NL) -Y2 (2,NL) ) * (Y  2 (1 ,NL) - 
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A2=Y2(1,NL)*Y2(3,NL)/((Y2(2,NL)-Y2(1,NL))*(Y2(2,NL)- 

Y  2  (  3  ,  NL  ) ) ) 

A3=Y2 ( 1 , NL ) *Y2 ( 2 , NL) / ( (Y2 ( 3 , NL ) -Y2 ( 1 , NL) ) * (Y2 ( 3 , NL)  - 

Y2  (2 ,NL) ) ) 


SUBSTITUTE  IN  INTERPOLATION  FORMULA 
DO  26  J= 1 , NL 

V1(J)-A1*Y2(1, J)+A2*Y2(2, J)+A3*Y2(3,J) 
V2(J)=A1*YY(1,J)+A2*YY(2, J ) +A3* YY ( 3 , J ) 

V3 (J)=A1*YYY(1 , J)+A2*YYY(2 , J)+A3*YYY (3 , J) 

P1(J)=A1*P(1,J) +A2*P (  2  ,  J  ) +A3*P (  3  ,  J  ) 

T1 ( J) =A1*T ( 1 , J) +A2*T (  2  ,  J) +A3*T ( 3 , J) 

THETA= FLOAT ( J-l) / 2  . 

WRITE (6,102)  THETA, VI(J),  V2(J),  V3(J),  Pl(J),  Tl(J) 
CONTINUE 
GO  TO  28 
CONTINUE 
DO  29  J= 1 , NL 
THETA=FLOAT (J-l) /2  . 

WRITE (6 ,102)  THETA,  Y2(LLL,J),  YY (LLL , J )  YYY (LLL , J ) , 
P (LLL , J) , 

IT  (LLL  ,  J) 

29  CONTINUE 
CONTINUE 

F0RMAT(1H1,//,55X, ’PLANE  FLOW’,//) 

FORMAT (// ,40X, 'ALPHA='  ,D8.2,10X, ' N= '  ,D10 . 2  ,  / ) 

FORMAT (10X,F4 . 1 , 6X , 5D  20 .8) 

FORMAT (/ ,10X, 'THETA' ,17X, 'M' ,19X, 'MM' , 1 7X , 'MMM  ,18X, 


26 


27 


28 

100 

101 

102 

108 


1 , 19X , ' T  '  ,/) 
STOP 
END 


'P' 
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SUBROUTINE  FCT (X, Y  ,  DERY ) 

DOUBLE  PRECISION  Y(3),  X,  DERY  (  3  )  ,  Cl,  C2,  TN 
COMMON  TN 
DERY (1) =Y (2) 

DERY ( 2 ) =Y  (3) 

Cl=4.0*Y(l) **2+Y (2) **2 
C2=2 .  *Y (2)* (4 , *Y (1) +Y ( 3 ) ) 

DERY ( 3 )  =  1 . / (C1+(TN-1 . ) *Y (2)**2) * ( (4  .*TN**2-12 . *TN+4 . ) 
1*Y(2)*C1-((4.*TN**2-6.*TN+2.)*Y(1)+(TN-1.)*Y(3))*C2- 
2(TN-1.)*(TN-3.)*Y(2)*C2**2/(4.*C1)-(TN-1.)*Y(2)*(4* 

3Y (2)**2+4*Y ( 1 )  *Y (3) +Y (3)**2)) 

RETURN 

END 
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SUBROUTINE  OUTP ( X , Y , DERY , IHLF , ND IM , PRHT ) 

DOUBLE  PRECISION  PRMT(5),  Y(3),  DERY (  3)  ,  AUX(8,3), 
1BETA,  ALPHA,  TN ,  Yl,YYl,  YYY1,  YYYYl ,  BETL,  BETH, 

2D  EL  ,  VI  ( 3 1 )  ,  V2(31),  V3(31),  Pl(31),  Tl(31),  Al,  A2, 
3A3,  PO,  Y2  ( 3 , 31 )  ,  YY(3,31),  YYY(3,31),  p(3,31),  T(3,31) 
COMMON  TN,BETA,Yl,PO ,Y2 , YY , YYY , P , T , LL , LLL 
IF  (DABS (X-LL*PRMT(3) *5)  .LT.0.1D-3)  GO  TO  20 
GO  TO  24 
20  LL=LL+1 

Y2  (LLL ,LL)=Y (1) 

YY (LLL,LL)=Y(2) 

YYY  (LLL  LL ) =Y  (  3 ) 

P(LLL,LL)=2**((1.-TN) / 2  .  ) * ( 2  . /TN* ( 1 . -TN ) *Y ( 1 ) - 1 . / 
2(2.*TN)*(Y(3)+Y(2)*(TN-l.)/2./(4.*Y(l)**2+Y(2)**2) 

3*  (8  ,*Y (1) *Y (2)  +  2 ,*Y (2)*Y(3) ) ) ) * ( 4 . *Y ( 1 ) ** 2+Y ( 2 ) ** 2 ) 

4**  (  (TN-1 . ) / 2 . )-PO 

T  (LLL ,LL)=Y (2)*(4.*Y(1) **2+Y (2)**2)** ( (TN-1 . )/2 .) 

2*2**  (  (1  .-TN) 12  .  ) 

24  IF (DABS (X-PRMT (2) )  .LT  .  0  .  ID-4)  GO  TO  12 
GO  TO  13 

12  Y 1=Y ( 1 ) 

13  CONTINUE 
RETURN 
END 


